
#########################################################
# Handling TCGA datasets
#########################################################

getTCGASampleNameFromFileName<-function(fileNames){
  sampleNames<-c()
  for(fileName in fileNames){
    
    tokens<-unlist(strsplit(fileName, '-', fixed=T))
    pId<-tokens[3]
    tissue<-tokens[4]
    if(length(grep(pattern='01A', tissue)) > 0){
      sampleNames<-c(sampleNames, paste(pId, sep='', 'T'))
    }
    else{
      sampleNames<-c(sampleNames, paste(pId, sep='', 'N'))
    }
  }
  return(sampleNames)
}

# RNASeq read counts
read.tcga.rnaseq.counts<-function(fnames){
  result<-data.frame()
  for(fname in fnames){
    sampleName<-regmatches(fname, regexpr('TCGA.*-[0-9][0-9]A', fname))
    aCountDF<-read.table(file=fname, header=T, row.names=1, stringsAsFactors=F)
    if(nrow(result) > 0){
      if(nrow(result) == nrow(aCountDF) && all(rownames(result) == rownames(aCountDF))){
        result<-cbind(result, aCountDF[,1])
      }
      else{
        print('Row name(s) does\'t equals')
      }
    }
    else{
      result<-data.frame(aCountDF[,1], row.names=rownames(aCountDF))
    }
    colnames(result)[ncol(result)]<-sampleName
  }
  return(result)
}

omit.zeroSums<-function(d, byrow=TRUE){
  if(byrow){
    keep<-rowSums(d) != 0
    return(d[keep,])
  }
  else{
    keep<-colSums(d) != 0
    return(d[,keep])    
  }
}


# clustering should be applied to intrinsic genes
subtypeCluster<-function(d, k=4, outDir='./output'){
  
  # gene median centered
  d<-sweep(d,1, apply(d,1,median,na.rm=T))

  # make consenesus clusters
  require(ConsensusClusterPlus)
  results<-ConsensusClusterPlus(d, maxK=k, reps=1000, pItem=0.8, pFeature=1, title=outDir,clusterAlg="hc",distance="pearson",seed=1262118388.71279,plot="pdf")

  return(results)
}

# GEO ExpressionSet
GEO.exprset<-function(gseAcc){
  require(GEOquery)
  require(Biobase)
  require(affy)
  require(gcrma)
# require(hgu133ahsentrezgcdf) #cdfname="HGU133A_HS_ENTREZG"
# require(hgu133ahsentrezgprobe)
# require(hgu133ahsentrezg.db)
  # Load GSE series
  gse<-getGEO(gseAcc, GSEMatrix = F)
  # Make sure that all of the GSMs are from the same platform
  platforms<-sapply(GSMList(gse), function(gsm) Meta(gsm)$platform)
  if(!all(platforms == platforms[1])){
    cat('Waring: there are different platforms from ', platforms[1])
  }
  else{
    # Get the probeset ordering
    probesets<-Table(GPLList(gse)[[1]])$ID    
  }
  # Make the data matrix from the VALUE columns from each GSM
  data.matrix<-do.call(cbind, lapply(GSMList(gse), function(x){
                                                    tab<-Table(x)
                                                    mymatch<-match(probesets, tab$ID_REF)
                                                    return(tab$VALUE[mymatch])
                                                }))
  data.matrix<-apply(data.matrix, 2, function(x) as.numeric(as.character(x)))
  data.matrix<-log2(data.matrix)
  # Go through the necessary steps to make a compliant ExpressionSet
  rownames(data.matrix)<-probesets
  colnames(data.matrix)<-names(GSMList(gse))
  pdata<-data.frame(samples=names(GSMList(gse)))
  rownames(pdata)<-names(GSMList(gse))
  pheno<-as(pdata, "AnnotatedDataFrame")
  return(new('ExpressionSet', exprs=data.matrix, phenoData=pheno))
}

mergePairedReadCounts<-function(subColNames, tumorReadCounts, normalReadCounts){
  result<-data.frame()
  for(subColName in subColNames){
    tumorColName<-paste(subColName, sep='', 'T')
    tumorColIndex<-match(tumorColName, colnames(tumorReadCounts))
    normalColName<-paste(subColName, sep='', 'N')
    normalColIndex<-match(normalColName, colnames(normalReadCounts))
    if(ncol(result) > 0){
      result<-data.frame(result, tumorReadCounts[, tumorColIndex], normalReadCounts[, normalColIndex], check.names=F)
    }
    else{
      result<-data.frame(tumorReadCounts[, tumorColIndex], normalReadCounts[, normalColIndex])
    }
    nc<-ncol(result)
    colnames(result)[c(nc-1, nc)]<-c(tumorColName, normalColName)
  }
  rownames(result)<-rownames(tumorReadCounts)
  return(result)
}

###############################################
## DEG analysis
###############################################
# At least one count per million (cpm) in all samples
filter.counts<-function(pairedCounts){
  nPatient<-ncol(pairedCounts)/2
  keep<-rowSums(cpm(pairedCounts[, grep('T', colnames(pairedCounts))]) > 1) >= nPatient | 
        rowSums(cpm(pairedCounts[, grep('N', colnames(pairedCounts))]) > 1) >= nPatient
  return(pairedCounts[keep,])  
}

# TMM normalization of read counts
tmm.normalize<-function(readCounts){
  require(edgeR)
  
  # create edgeR object and TMM normalization
  d<-DGEList(targetReadCounts)
  d<-calcNormFactors(d, method='TMM')
  
  # calculate log-counts-per-million
  return(cpm(d, normalized.lib.sizes=TRUE, log=TRUE))
}

# Estimate DEGs using edgeR
# Create negative binomial GLM and perform LRT
# Column names of input data are SampleId + ('T'| 'N'), eg., "6771T" "6771N" "5472T" "5472N"...
fit.glm<-function(pairedCounts){
  # Create edgeR object and normalization
  y<-DGEList(counts=pairedCounts)
  y<-calcNormFactors(y, method='TMM')

  # Experimental design matrix
  cns<-colnames(y)
  Patient<-factor(substring(cns, 1, 4))
  Tissue<-factor(substring(cns, 5, 5))
  design<-model.matrix(~Patient+Tissue)
  rownames(design)<-cns
  
  ## Estimate dispersion
  y <- estimateGLMCommonDisp(y, design, verbose=TRUE) # Estimates common dispersions
  y <- estimateGLMTrendedDisp(y, design) # Estimates trended dispersions
  y <- estimateGLMTagwiseDisp(y, design) # Estimates tagwise dispersions 
  
  ## Fit the negative binomial GLM for each tag
  fit <- glmFit(y, design) # Returns an object of class DGEGLM
  lrt <- glmLRT(fit)
  
  return(lrt)
}

######################################################
## Network analysis
######################################################
## Build differential expressed gene weighted graph
build.degGraph<-function(eTab, vTab){
  require(igraph)
  g<-graph.data.frame(eTab, directed=F, vertices=vTab)
  
  # calculate edge weights using vertex's logFC
  E(g)$Weight<-0
  for(v in V(g)){
    #cat('===========>Vertex:', V(g)[v]$name, '\n')
    cat('===========>Vertex:', v, '\n')
    
    incidents<-E(g)[inc(v)]
    if(length(incidents) > 0){
      vw<-vertex.weight(g, v)
      cat('Vertex weight:', vw, '\n')
      E(g)[incidents]$Weight<-E(g)[incidents]$Weight + (vw/length(incidents))
      #cat('Incident edges:', incidents, '\n')
      #cat('Incident weights:', E(g)[incidents]$Weight, '\n')
    }    
  }
  return(g)
}

vertex.weight<-function(g, v){
  return(abs(V(g)[v]$logFC))
}

get.moduleIndices<-function(vTab, modules){
  result<-vTab
  moduleIndices<-rep(NA, nrow(vTab))
  for(i in seq(1:length(modules))){
    found<-match(modules[[i]], rownames(vTab))
    moduleIndices[found]<-paste2(moduleIndices[found], i, sep=',')
  }
  return(moduleIndices)
}

build.subGraphs<-function(vTab, eTab, modules){
  require(igraph)
  results<-list()
  for(module in modules){
    subVTab<-vTab[vTab$name %in% module,]
    subETab<-eTab[(eTab$from %in% module)&(eTab$to %in% module),]
    subGraph<-graph.data.frame(subETab, directed=F, vertices=subVTab)
    subGraph$PertScore<-mean(subETab$Weight)
    results[[length(results) + 1]]<-subGraph
  }
  
  return(results)
}

################################################################
# Pathway Analysis
################################################################
# Fisher's exact test
# getContingencyMatrix(module, length(ppi.genes$ID), length(deg$GeneId))
getContingencyMatrix<-function(geneIds, N, n){
  r<-length(geneIds)
  a<-length(which(!is.na(match(geneIds, oscc.filteredDETable$GeneId))))
  b<-n-a
  c<-r-a
  d<-N-r-b
  if(a > 0){
    mm<-matrix(c(a, b, c, d), nrow=2, dimnames=list(c('InModule', 'NotInModule'),c('DEG', 'NotDEG')))
    mm
  }
}

reactome.pathway.enrich<-function(pathways, geneSets, cutoff=0.05, N=30000){
  results<-list()
  for(genes in geneSets){
    qvals<-qval.pathway.enrich(pathways, genes, N)
    item<-list(Genes=genes, ReactomeIDs=sapply(pathways[qvals<=cutoff], function(pw) pw$ReactomeID))
    results[[length(results) + 1L]]<-item
  }
  return(results)
}

## Get qvalues(adjusted pvalues) to test whether each pathway significantly contains the genes(symbols)
# using Fisher's exact test
## Args
# pathways  a list of pathways, each pathway has 'Name', 'Genes' attributes
# symbols a vector of gene symbol to test significancy
# N the human genome background (30,000 gene total)
## Values
# a vector of q values which are adjusted p values by FDR
qval.pathway.enrich<-function(pathways, genes, N=30000){
  pvals<-sapply(pathways, function(pathway) pval.overlap.genesets(pathway$Genes, genes, N))
  
  # correct multiple comparison problem, adjust p-values
  return(p.adjust(pvals, method='fdr'))
}


# Test whether source genes significantly overlap with target genes,
# Using Fisher's exact test
# N: Total number of genes in human genome
pval.overlap.genesets<-function(targetGenes, sourceGenes, N=30000){
  
  nt<-length(sourceGenes)
  np<-length(targetGenes)
  
  a<-length(which(sourceGenes %in% targetGenes))
  b<-nt - a
  c<-np - a
  d<-N - np - b
  
  m<-matrix(c(a, b, c, d), nrow=2, dimnames = list( c('InSource', 'NotInSource'), c('InTarget', 'NotInTarget')), byrow=T)
  
  ft<-fisher.test(m, alternative = 'greater')
  
  return(ft$p.value)	
}
#################################################################
# Handling CTRP datasets
# Ref: Basu, A. et. al., An Interactive Resource to Identify Cancer Genetic and Lineage Dependencies Targeted by Small Molecules. Cell, 2013.
# Suppl. DataSet: http://www.cell.com/cms/attachment/2009208227/2031783034/mmc2.xlsx
# Note: The suppl. data file(mmc2.xlsx) is devid into sub files, e.g., sub_s2b.xlsx, sub_s2g.xlsx
# whose names are given by the Suppl. data user guide(http://www.cell.com/cms/attachment/2009208227/2031783072/mmc6.pdf)
#################################################################
read.ds<-function(){
  require(xlsx)
  df<-read.xlsx('datasets//ctrp/sub_s2d.xlsx', sheetIndex = 1, header=T, check.names=F)
  colnames(df)[1]<-'GeneName'
  df$GeneName<-as.character(df$GeneName)
  tmp<-geneID.from(df$GeneName)
  # df$GeneName[665] is 'LOC100130776', whose ID is not given by geneID.from
  tmp[665]<-'100130776'
  rownames(df)<-tmp
}

#################################################################
# Common Utility Functions
#################################################################
substr2<-function(x, sep=',', idx){
  return(sapply(x, function(cur) unlist(strsplit(cur, sep, fixed = T))[idx]))
}
paste2 <- function(..., sep='') {
  L <- list(...)
  L <- lapply(L,function(x) {x[is.na(x)] <- ""; x})
  ret <-gsub(paste0("(^",sep,"|",sep,"$)"),"",
             gsub(paste0(sep,sep),sep,
                  do.call(paste,c(L,list(sep=sep)))))
  is.na(ret) <- ret==""
  ret
}

geneID.from<-function(symbols){
  require(org.Hs.eg.db)
  idTable<-toTable(org.Hs.egSYMBOL2EG)
  
  idx<-match(toupper(symbols), toupper(idTable$symbol))
  found<-idTable$gene_id[idx]
  
  idTable<-toTable(org.Hs.egALIAS2EG)
  
  idx<-match(toupper(symbols[is.na(found)]), toupper(idTable$alias_symbol))
  
  found[is.na(found)]<-idTable$gene_id[idx]
  return(found)
  
}  

read.reactome.pathwayGenes<-function(fname){
  lines<-readLines(fname)
  return(lapply(strsplit(lines, '\t', fixed=T), function(x) list(Name=x[1], Genes=x[3:length(x)])))
}

read.list<-function(fname, cutoff=1, callback.lineItems=function(x) x){
  results<-list()
  for(line in readLines(fname)){
    items<-unlist(strsplit(line, '\t', fixed=T))
    if(length(items) > cutoff){
      results<-append(results, list(Items=callback.lineItems(items)))
    }
    
  }
  return(results)
}

convertReactomeSymbolsIfPossible<-function(symbols){
  	found<-match(symbols, reactome.fi.genes$ReactomeSymbol)
  	target<-reactome.fi.genes$EntrezSymbol[found]
	target[is.na(found)]<-symbols[is.na(found)]
  	return(target)
}

reactomeSymbol2EntrezID<-function(symbols){
  return(unique(na.omit(reactome.fi.genes$EntrezID[match(symbols, reactome.fi.genes$ReactomeSymbol)])))
}

